import os
import sys
import shutil
import logging
import glob
import matplotlib.pyplot as plt
from matplotlib import rc
rc('text', usetex = True)
plt.rcParams.update( # try to match font sizes of document
{'axes.labelsize': 20,
'axes.titlesize': 20,
'legend.fontsize': 20,
'xtick.labelsize': 20,
'ytick.labelsize': 20,
'text.usetex': True,
'font.family': 'serif',
'font.serif': ['palatino'],
'savefig.dpi': 300
})
%pylab inline
from glue.ligolw import ligolw
from glue.ligolw import table
from glue.ligolw import lsctables
from glue.ligolw import ilwd
from glue.ligolw import utils as ligolw_utils
import pycbc.strain
import pycbc.psd
from pycbc.pnutils import mass1_mass2_to_mchirp_eta
from pycbc.waveform import td_approximants, fd_approximants
from pycbc.waveform import get_two_pol_waveform_filter, get_td_waveform
from pycbc import DYN_RANGE_FAC
from pycbc.types import FrequencySeries, zeros
from pycbc.filter import match, overlap, sigma, make_frequency_series
from pycbc.scheme import CPUScheme, CUDAScheme
from GWNRTools.Utils.SupportFunctions import make_padded_frequency_series
from pycbc import pnutils
sys.path.append('/home/prayush.kumar/local/venv/pycbc_master_enigma/src/GWNRTools/bin/')
run_dirs = [
'/home/prayush.kumar/projects/template_banks/faithfulness_SEOBNRv4_vs_SEOBNRv4_ROM',
'/home/prayush.kumar/projects/template_banks/faithfulness_SEOBNRv4_vs_SEOBNRv2_ROM',
'/home/prayush.kumar/projects/template_banks/faithfulness_SEOBNRv4_vs_SEOBNRv2',
'/home/prayush.kumar/projects/template_banks/faithfulness_SEOBNRv4_vs_EOBNRv2'
]
models = [
'SEOBNRv4_ROM', 'SEOBNRv2_ROM_DoubleSpin', 'SEOBNRv2', 'EOBNRv2'
]
plot_dir = '/home/prayush.kumar/src/EccentricIMRTemplateBank/Paper/plots/faithfulness/'
for d in run_dirs:
try: os.makedirs(d)
except: logging.warn("Directory {} already exists.".format(d))
os.chdir()
for d in run_dirs:
os.chdir(d)
!ls
!rm -fr *
!rm -rf bank log match plots scripts *.sub faithsim.dag* faithsim.sh inj*xml dag.out result*.dat submit.sh
!ls
ini_filename = "faithsim.ini"
for d, approx in zip(run_dirs, models):
os.chdir(d)
with open(ini_filename, "w") as fp:
fp.write("""
[inspinj]
num-new-points = 5000
component-mass-min = 5.0
component-mass-max = 50.0
total-mass-max = 100.0
spin-component-min = 0
spin-component-max = 0
spin-mag-min = 0
spin-mag-max = 0
eccentricity-min = 0
eccentricity-max = 0
coa-phase-min = 0
coa-phase-max = 0
inclination-min = 0
inclination-max = 0
long-asc-nodes-min = 0
long-asc-nodes-max = 0
mean-per-ano-min = 0
mean-per-ano-max = 0
mchirp-window = 0
eccentricity-window = 0
output-prefix = inj
verbose =
[executables]
inspinj = /home/prayush.kumar/local/venv/pycbc_master_enigma/src/GWNRTools/bin/gwnrtools_sample_parameter_space
faithsim = /home/prayush.kumar/local/venv/pycbc_master_enigma/src/GWNRTools/bin/gwnrtools_faithsim
[workflow]
accounting-group = ligo.dev.o3.cbc.explore.test
templates-per-job = 100
log-path = /usr1/prayush.kumar/
[faithsim-flatIMRC]
psd-model = aLIGOZeroDetHighPower
waveform1-approximant = {0}
waveform1-start-frequency=24
waveform2-approximant = SEOBNRv4
waveform2-start-frequency=24
filter-low-frequency=25
filter-sample-rate=8192
filter-waveform-length=64
""".format(approx))
for d in run_dirs:
os.chdir(d)
logging.warn("Configuring faithsim run in {0}".format(d))
!/home/prayush.kumar/local/venv/pycbc_master_enigma/src/GWNRTools/bin/gwnrtools_create_faithsim_workflow \
--config=faithsim.ini
def param_from_table(param, table):
return np.array([
getattr(p, param) for p in table
])
for d, approx in zip(run_dirs, models):
inj_filename = 'inj.xml'
indoc = ligolw_utils.load_filename(inj_filename, contenthandler=table.use_in(ligolw.LIGOLWContentHandler))
inj_table = lsctables.SimInspiralTable.get_table(indoc)
print("Total {0} test-points in {1}".format(len(inj_table), d))
m1 = np.array([inj.mass1 for inj in inj_table])
m2 = np.array([inj.mass2 for inj in inj_table])
plt.figure()
plt.plot(m1, m2, 'kx', alpha=0.3)
plt.grid()
plt.title(approx.replace('_', '-'))
plt.xlabel('mass 1')
plt.ylabel('mass 2')
for d in run_dirs:
!ls
for d in run_dirs:
os.chdir(d)
logging.warn("Submitting workflow in {0}".format(d))
!condor_submit_dag faithsim.dag >> dag.out
!cat dag.out
for d in run_dirs:
os.chdir(d)
logging.warn("Status of workflow in {0}".format(d))
!tail faithsim.dag.dagman.out
bfields = ('match', 'overlap', 'time_offset', 'sigma1', 'sigma2', 'mass1',
'mass2', 'spin1x', 'spin1y', 'spin1z', 'spin2x', 'spin2y',
'spin2z', 'inclination', 'latitude', 'longitude',
'polarization', 'coa_phase')
dtypeo={'names': bfields,
'formats': ('f8', 'f8', 'f8', 'f8', 'f8', 'f8', 'f8', 'f8', 'f8',
'f8', 'f8', 'f8', 'f8', 'f8', 'f8', 'f8', 'f8', 'f8')}
all_results = {}
all_masked_results = {}
for d, approx in zip(run_dirs, models):
os.chdir(d)
results_filename = glob.glob('result-*.dat')[0] #'result-flatIMRC.dat'
all_results[approx] = np.loadtxt(results_filename, dtype = dtypeo)
# clean up data results
all_masked_results[approx] = all_results[approx][np.isfinite(all_results[approx]['match'])]
all_masked_results[approx] = all_masked_results[approx][all_masked_results[approx]['match'] > 0.0]
print(len(all_masked_results[approx]['match']))
print(all_masked_results[approx]['match'])
for d, a in zip(run_dirs, models):
logging.warn("Showing figures for {0}".format(a))
masked_results = all_masked_results[a]
# figure style 1
plt.figure(figsize = (8, 5))
plt.scatter(masked_results['mass1'], masked_results['mass2'], c = masked_results['match'],
s = 20, vmax = 1.0)#, vmin = 0.96)
plt.xlabel('$m_1 (M_\odot)$', fontsize = 20)
plt.ylabel('$m_2 (M_\odot)$', fontsize = 20)
plt.grid()
cb = plt.colorbar()#(label='Overlap')
cb.set_label('Match vs {}'.format(a.replace('_', '\_')))
plt.savefig(os.path.join(plot_dir,
'{0}_SEOBNRv4_no_incl_no_ecc_match_over_mass1_mass2_scatter.pdf'.format(a)))
# figure style 2
plt.figure(figsize = (8, 5))
plt.hexbin(masked_results['mass1'], masked_results['mass2'], C = masked_results['match'],
gridsize = 40, vmax = 1.0)
plt.xlabel('$m_1 (M_\odot)$', fontsize = 20)
plt.ylabel('$m_2 (M_\odot)$', fontsize = 20)
plt.grid()
cb = plt.colorbar()#(label='Overlap')
cb.set_label('Match vs {}'.format(a.replace('_', '\_')))
plt.savefig(os.path.join(plot_dir,
'{0}_SEOBNRv4_no_incl_no_ecc_match_over_mass1_mass2_hex.pdf'.format(a)))
for d, a in zip(run_dirs, models):
logging.warn("Showing figures for {0}".format(a))
masked_results = all_masked_results[a]
# figure style 1
plt.figure(figsize = (8, 5))
plt.scatter(masked_results['mass1'], masked_results['mass2'], c = 1. - masked_results['match'],
s = 20, vmin = 0)#, vmin = 0.96)
plt.xlabel('$m_1 (M_\odot)$', fontsize = 20)
plt.ylabel('$m_2 (M_\odot)$', fontsize = 20)
plt.grid()
cb = plt.colorbar()#(label='Overlap')
cb.set_label('Mismatch vs {}'.format(a.replace('_', '\_')))
plt.savefig(os.path.join(plot_dir,
'{0}_SEOBNRv4_no_incl_no_ecc_mismatch_over_mass1_mass2_scatter.pdf'.format(a)))
# figure style 2
plt.figure(figsize = (8, 5))
plt.hexbin(masked_results['mass1'], masked_results['mass2'], C = 1. - masked_results['match'],
gridsize = 40, vmin = 0)
plt.xlabel('$m_1 (M_\odot)$', fontsize = 20)
plt.ylabel('$m_2 (M_\odot)$', fontsize = 20)
plt.grid()
cb = plt.colorbar()#(label='Overlap')
cb.set_label('Mismatch vs {}'.format(a.replace('_', '\_')))
plt.savefig(os.path.join(plot_dir,
'{0}_SEOBNRv4_no_incl_no_ecc_mismatch_over_mass1_mass2_hex.pdf'.format(a)))
# figure style 3
plt.figure(figsize = (8, 5))
plt.hexbin(masked_results['mass1'], masked_results['mass2'], C = 1. - masked_results['match'],
gridsize = 40, vmax = 0.01, vmin = 0.)
plt.xlabel('$m_1 (M_\odot)$', fontsize = 20)
plt.ylabel('$m_2 (M_\odot)$', fontsize = 20)
plt.grid()
cb = plt.colorbar()#(label='Overlap')
cb.set_label('Mismatch vs {}'.format(a.replace('_', '\_')))
#plt.savefig(os.path.join(plot_dir, 'SEOBNRv4_ENIGMA_no_incl_no_ecc_match_over_mass1_mass2_hex.pdf'))
for d, a in zip(run_dirs, models):
logging.warn("Showing figures for {0}".format(a))
masked_results = all_masked_results[a]
# figure style 1
plt.figure(figsize = (8, 5))
plt.scatter(masked_results['mass1'], masked_results['mass2'], c = masked_results['match'],
s = 20, vmax = 1)#, vmin = 0.96)
plt.xlabel('$m_1 (M_\odot)$', fontsize = 20)
plt.ylabel('$m_2 (M_\odot)$', fontsize = 20)
plt.grid()
cb = plt.colorbar()#(label='Overlap')
cb.set_label('Match vs {}'.format(a.replace('_', '\_')))
plt.savefig(os.path.join(plot_dir,
'{0}_SEOBNRv4_no_incl_no_ecc_match_over_mass1_mass2_scatter.pdf'.format(a)))
# figure style 2
plt.figure(figsize = (8, 5))
plt.hexbin(masked_results['mass1'], masked_results['mass2'], C = masked_results['match'],
gridsize = 40, vmax = 1)
plt.xlabel('$m_1 (M_\odot)$', fontsize = 20)
plt.ylabel('$m_2 (M_\odot)$', fontsize = 20)
plt.grid()
cb = plt.colorbar()#(label='Overlap')
cb.set_label('Match vs {}'.format(a.replace('_', '\_')))
plt.savefig(os.path.join(plot_dir,
'{0}_SEOBNRv4_no_incl_no_ecc_match_over_mass1_mass2_hex.pdf'.format(a)))
# figure style 3
plt.figure(figsize = (8, 5))
plt.hexbin(masked_results['mass1'], masked_results['mass2'], C = masked_results['match'],
gridsize = 40, vmax = 1, vmin = 0.99)
plt.xlabel('$m_1 (M_\odot)$', fontsize = 20)
plt.ylabel('$m_2 (M_\odot)$', fontsize = 20)
plt.grid()
cb = plt.colorbar()#(label='Overlap')
cb.set_label('Match vs {}'.format(a.replace('_', '\_')))
#plt.savefig(os.path.join(plot_dir, 'SEOBNRv4_ENIGMA_no_incl_no_ecc_match_over_mass1_mass2_hex.pdf'))